Efficient Monte Carlo Calculations of the One-Body Density 
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Abstract 

An alternative Monte Carlo estimator for the one-body density p(r) is presented. This estimator 
has a simple form and can be readily used in any type of Monte Carlo simulation. Comparisons 
with the usual regularization of the delta-function on a grid show that the statistical errors are 
greatly reduced. Furthermore, our expression allows accurate calculations of the density at any 
point in space, even in the regions never visited during the Monte Carlo simulation. The method 
is illustrated with the computation of accurate Variational Monte Carlo electronic densities for 
the Helium atom (ID curve) and for the water dimer (3D grid containing up to 51x51x51=132651 
points). 

PACS numbers: 02.70.Ss, 02.70.Uu, 02.50.Ng,71.15.-m 
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The Monte Carlo approach is probably one of the most widely employed numerical ap- 
proaches in the scientific and engineering community. In computational physics, it has been 
extensively used in the last fifty years for studying a great variety of many-body systems un- 
der many different conditions. To date, the most popular application of the method is prob- 
ably the calculation of classical thermodynamical properties, pj However, the Monte Carlo 
approach is also employed for evaluating quantum properties by using the Path-Integral 
formulation of quantum averages as classical ones (Quantum Monte Carlo or Path Integral 
Monte Carlo approaches [2]). In the recent years, these later approaches have emerged as 
an unique and powerful tool for studying quantitatively the interplay between quantum and 
thermal effects in many-body systems (e.g., to understand the very rich physics of strongly 
correlated materials). 

At the heart of all these applications lies the calculation of a number of high-dimensional 
integrals (or sums, for lattice problems) written under the general form 

1(F) = j dr x ■ ■ ■ dr N U( ri , ■ ■ • , r N )F(n, ■ ■ ■ , r N ) (1) 

where II is some arbitrary iV-body probability distribution (II positive and normalized) and 
F some arbitrary real-valued function. The integration is performed over all accessible con- 
figurations for the iV-particle system. The general idea of Monte Carlo approaches is to 
evaluate the integral by sampling the configuration space according to the probability dis- 
tribution, II, and by averaging F over the various configurations generated by the sampling 
procedure, 1(F) = (F)jj. Here and in what follows, the symbol (• • -)n indicates the statistical 
average over the density n. Various Monte Carlo algorithms (sampling procedures) can be 
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found in the literature, the most celebrated one being, of course, the Metropolis algorithm 
The efficiency of a Monte Carlo approach is directly related to the magnitude of the fluctu- 
ations of the integrand in the regions where the probability distribution, n, is large. More 
precisely, for a given number of Monte Carlo steps, the statistical error SF is proportional 
to the square root of the variance of the integrand F defined as <J 2 (F) = I(F 2 ) — 1(F) 2 . 
Accordingly, a very attractive way of enhancing the convergence of a Monte Carlo simula- 
tion consists in introducing alternative "improved" estimators defined as new integrands F 
having the same average as F but a lower variance 

(F)n=(F)n and a 2 (F)<a 2 (F). (2) 
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In previous worksP, it has been shown how improved estimators can be designed for any 
type of integrand F and Monte Carlo algorithm, and some applications to the computation 
of forces have been presented . 

In this Letter we present an efficient improved Monte Carlo estimator for calculating the 
one-body (or one-particle) density, p(r) 

N 

p(r) = <£*(ri-r)>n (3) 

i=l 

and, more generally, any one-body average of the form J drp(r)F(r). As we shall see, 
our estimator allows very important reductions in variance. In the example of the charge 
density of the water dimer presented below, a reduction of up to two orders of magnitude 
in CPU time is possible for some regions of space. In addition, and in sharp constrast 
with the usual estimator based on the regularization of the delta- function on a grid, our 
expression leads to accurate estimates of the density at any point in space, even in the regions 
never visited during the Monte Carlo simulation (e.g., in the large-distance regime). This 
property is particularly interesting when a global knowledge of the density map is searched 
for. For the water dimer case, we were able to accurately compute the charge density for 
51x51x51=132651 grid points. Note that such a calculation is vastly more difficult to perform 
when using the standard approach. 

Let us recall that accurate one-particle properties are of central interest for the under- 
standing of the physics of many complex many-body systems. Such systems include all 
those which are not translationally invariant (typically, all finite systems: atoms, molecules, 
clusters, nuclei, etc.) and all those whose translational symmetry has been explicitly broken, 
e.g., by the application of an inhomogeneous external field. In addition to this, many phys- 
ical modelizations and/or effective theories rely explicitly on the knowledge of the one-body 
density. Many examples could be cited but let us mention, for example, the various mod- 
elings of the electrostatic field of molecules from the one-electron density 3], the studies 
of the structure and reactivity of molecular systems based on the topological analysis of 
the electron density and/or its Laplacian, Q| , and, also, the very important case of Density 
Functional Theories (DFT) which could greatly benefit from the possibility of computing 
accurate 3D charge/spin density maps for large molecular systems (e.g., via accurate fits of 
the exchange-correlation Kohn-Sham potential, see ^^). 

Finally, let us note that the use of alternative forms for evaluating the density is not 



new. For example, in the works of Hiller et aZ.[f|, Sucher and DrachmanQ, Harimanjsj], 
Rassolov and Chipman|9] new classes of global operators built for computing the density 
have been introduced. However, in these works, the general idea is to design operators 
whose expectation values give an accurate estimate of the unknown exact one-body density 
and not the exact one-body density associated with a known A-body density. Actually, our 



strategy is more closely related to w 



rat has been presented by Vrbik et al. |l0( , Langfelder et 



a/.[n|, and Alexander and Coldwell 3]. In these works, alternative Monte Carlo estimators 
with lower variances are also introduced. However, the emphasis is only put on the case of 
evaluating the charge and/or spin density at the nuclei. Here, such ideas are extended to 
any point in space and a general formula allowing to control all possible sources of statistical 
fluctuations in all possible regimes is presented. 

General improved density estimator. Due to the presence of Dirac functions in the Monte 
Carlo estimator of the density, Eq.Q, some sort of regularization has to be introduced. It is 
usually done by partionning the physically relevant part of the one-particle space (usually, 
the 3D ordinary space) into small domains of finite volume and by evaluating the corre- 
sponding locally-averaged densities. In practice, such a procedure is particularly simple to 
implement by counting the number of particles present in each elementary domain at each 
step of the simulation. However, the statistical fluctuations can be rather large. This is 
particularly true for the low-density regions which are rarely visited by the particles. Even 
worse, there is no way of evaluating the density in regions which are never visited during the 
finite Monte Carlo simulation. One way to escape from these difficulties is to introduce some 
global estimators defined in the whole space. To do that, we regularize the 5 Dirac-function 
by using the following equality 

^-H^V^, (4) 

where /(ri; r) is a smooth function of ri (here, r plays the role of an external parameter) 
verifying /(r$ = r; r) = 1. Note that this formula is just a slightly generalized form of the 
well-known equality corresponding to / = 1. Now, injecting this expression into Eq.(J3J) and 
integrating by parts, the density can be rewritten as 

i=i 1 * 1 

Next, we introduce some additional function g(r) independent on the particle coordinates 



and write p(r) under the form 

P r = I] 1 1 - 9} n n, 6 

47T Tj — r 11 

i=l 1 1 

v 2 f fro 

the last step being allowed since ( ^ )n = 0. 

Expression (JHJ) is our general form for the improved estimator of the one-body density. 
The two functions / and g play the role of auxiliary quantities. They are introduced to 
decrease the variance of the density estimator. As with any optimization problem, there is 
no universal strategy for choosing / and g. However, the guiding principle is to identify the 
leading sources of fluctuations and, then, to adjust the auxiliary functions to remove most 
of them. 

Improved electronic density estimator for molecules. 

In what follows, we consider the one-electron density of molecular systems issued from 
the iV-body quantum probability distribution written as 

n = iPU r u---,r N ) ^ 

/ dri - ■ -drx-ip^T!, - ■ ■ ,r N )' 

where ipT is some electronic trial wavefunction. 
1. Short electron-nucleus distance regime. 

For a system of electrons in coulombic interaction with a set of fixed nuclei the exact 
wave function is known to obey the following electron-nuclear cusp condition 

V> ~r,^R A 1 - Z A \vi - Ra\, (8) 

where R4 denotes the position of a given nucleus A of charge Z A . Most of the accurate 
trial wavefunctions employed in the literature fulfill this important condition. Now, as a 
consequence of the cusp condition we have 

v 2 n az a 
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(9) 



In the neigbhorhood of nucleus A, this term is an important source of fluctuations. This is 

easily seen by noting that, in the regime r ~ R^, the estimator of p(r), Eq.flOJ), behaves as 

|r--R A | 2 ' a Q uan tity which has an infinite variance [J (^-) 2 r 2 (ir = +00]. To remove this source 

v 2 f 

of wild fluctuations we adjust the function / so that — y- exactly vanishes the divergence 
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of A simple suitable form for / verifying such a condition plus the constraint, /(r, = 
r; r) = 1, is given by 

/(r i; r) = 1 + 2Z A (\ ri - R A \ - |r - R A |). (10) 
2. Large-distance regime. 

In the large-distance regime, |r| — > +oo, the exact one-electron density is known to decay 
exponentially. In contrast, our basic estimator decays algebraically as a function of r. Here, 
we propose to force the latter estimator to decay also exponentially. In practice, a simple 
choice for the function / is 

/(r i; r) = (1 + A|r, - r|) exp [-A|r, - r|], (11) 

where A is some real parameter. Note that the coefficients of the linear and exponential terms 
have been taken identical (= A) to avoid the divergence of the laplacian of / at = r. The 
parameter A can be adjusted eiher by minimizing the fluctuations of the average density 
or fixed at some value close to the theoretical value of X ex = 2y '—21 (I first ionization 
potential, the exact density is known to decay as p ~ exp — X ex r). 

Another useful remark is that, in the regime |r| — > +oo, the electrons of the molecule 
are, in first approximation, at the same distance from the point where p is evaluated. As 
a consequence, during the simulation the quantity 1/ — r| fluctuates very little around 
the approximate average 1/j (r») — r|. Therefore, to reduce the fluctuations it is valuable 
to remove a quantity close to the latter average from the former one. Here, this idea is 
implemented by introducing the following function g 

M 

A=l 1 A 1 

In our first application we consider the He atom described by a simple trial wavefunction 
written under the form, ip T = 0(r x )0(r 2 ) with = exp (— jr), and 7 = 1.6875 (Slater value). 
For this problem, the exact density is known and is given by p(r) = 2 / y 3 / tc exp (—277*). Figure 
n shows the results obtained for p(r) for a relatively short Monte Carlo run. The main curve 
displays the results obtained with i.) the usual estimator based on the delta representation, 
Eq.Q, ii.) the simple improved estimator corresponding to Eq.© with / = 1 and g = 0, 
and our best improved estimator defined via Eqs. (j6ll0lllll2J) . In the latter case, the densities 
corresponding to each of the two possible choices for / [Eqs.(JTUJ) or (JTTJ)] have been computed 
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for each distance, the final value corresponding to the value having the smallest statistical 
error. The usual estimator, Eq.(j3J), has been regularized by introducing small elementary 
cubes of length a = 0.2. For all distances the statistical error associated with the usual 
estimator is very large with respect to improved estimators. At intermediate distances, at 
least one order of magnitude in accuracy is lost. For example, at r = 0.6 the statistical error 
is about 10 times larger than for the simple estimator case and a factor of about 20 is found 
with respect to the best improved estimator. At large distances, a region rarely visited by 
the electrons, the standard estimator is so noisy that it is useless in practice. Now, regarding 
improved estimators it is clear that the auxiliary functions / and g have a great impact in 
reducing the errors. At very small distances (first inset) a gain of about 5 in statistical error 
is obtained with the best improved estimator. At r = 0, the gain is even larger since the 
simple estimator has an infinite variance. At large distances where both improved estimators 
have a finite variance, it is seen that introducing some exponential decay into the estimator 
plus a proper shift ((^-contribution) improves considerably the convergence. In the range 
(2.5,3.) the gain in error increases from 15 at r = 2.5 to 40 at r = 3. 

In our second application we consider the water dimer in a non-symmetric nuclear ge- 
ometry (structure #2 of Ref. Q|) described by an electronic wavefunction consisting of a 
Hartree-Fock part (cc-pVTZ basis set) plus a standard explicitly correlated Jastrow term. 
Figure El shows the density plots obtained with our best improved estimator for a number 
of points equal to 51x51x51 (=132651). As seen on the figure the density obtained displays 
a very smooth aspect. A closer look shows that this regularity is present at a rather small 
scale. In Figure El we present a more quantitative comparison of the data along the 0-0 
axis. The figure clearly shows that the best estimator outperforms the usual one. First, the 
curve corresponding to the new estimator (solid line connecting the points) is very smooth, 
although it has been obtained by simple linear extrapolation of the data. In sharp con- 
trast, this is absolutely not true for the usual estimator curve whose overall behavior is 
particularly chaotic. Second, the statistical error has been greatly reduced using the new 
estimator. Depending on the distance, a gain in accuracy ranging roughly from 5 to 10 (i.e, 
up to two orders of magnitude in CPU time) has been obtained. An interesting point to 
mention is the presence of some very wild fluctuations in the neighborhood of r~ 1.5 for the 
standard estimator. These fluctuations are due to the presence of a hydrogen atom close 
to the — axis. We can verify that, in sharp contrast, our new estimator, which has 
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been built to correctly take into account the nuclear cusp, Eq.(|10|). performs well in that 
region. Finally, remark that in the large-r regime (data not shown here) where the standard 
estimator is strictly zero (no sampling of this region), the improved estimator still continues 
to give accurate values of the very small density 
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FIGURE CAPTIONS 



• FigC] Density of Helium from various estimators, see text. 

• FigJ2] One-electron density of the H 2 dimer with the best improved estimator. 

• FigEl Cut of the one-electron density along the 0-0 axis of the H 2 dimer. Data for 
the best and usual estimators. Solid lines are simple linear extrapolations of the data. 
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